(-1.5, 1.5)
(-1.5, 1.5)
Ingeniería Biomédica
2025-12-02
(Continuous-Time)
. . .
Signal Domain: Continuous signal: \(x(t)\)
. . .
Definition (Analysis): \[X(j\Omega) = \int_{-\infty}^{\infty} x(t) e^{-j\Omega t} dt\]
. . .
Frequency Variable: \(\Omega\) (Continuous frequency, [rad/s])
. . .
Spectrum Nature: The spectrum \(X(j\Omega)\) is: * Continuous * Aperiodic
. . .
. . .
Signal Domain: Discrete signal: \(x[n]\)
. . .
Definition (Analysis): \[X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}\]
. . .
Frequency Variable: \(\omega\) (Normalized frequency, [rad/sample])
. . .
Spectrum Nature: The spectrum \(X(e^{j\omega})\) is: * Continuous * Periodic (with period \(2\pi\))
. . .
The Bridge: The DTFT is the Z-Transform, \(X(z)\), evaluated precisely on the unit circle (\(z = e^{j\omega}\)).
Let \(x[n]\) be a discrete-time signal.
The Z-Transform is defined as:
\[X(z) = \sum_{n=-\infty}^{\infty} x[n] z^{-n}\]
Where: - \(z \in \mathbb{C}\) is a complex variable - \(z = re^{j\omega}\)
Causal Signals ROC is outside outermost pole
Anti-Causal Signals ROC is inside innermost pole
\[H(z) = 1.00 \cdot \frac{(z - 0.50)}{(z - 0.90)}\]
(-1.5, 1.5)
(-1.5, 1.5)
If the ROC includes the unit circle, \(|z| = 1\), then:
\[X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}\]
So the Fourier Transform is a special case of the Z-Transform.
Let:
\[x[n] = a^n u[n]\]
Where: - \(a \in \mathbb{R}\) - \(u[n]\) is the unit step function (0 for \(n<0\), 1 for \(n\geq 0\))
Apply the definition:
\[X(z) = \sum_{n=0}^{\infty} a^n z^{-n} = \sum_{n=0}^{\infty} (az^{-1})^n\]
This is a geometric series:
\[X(z) = \frac{1}{1 - az^{-1}} = \frac{z}{z - a}\]
A geometric series is a sum of terms where each term is multiplied by the same constant:
\[S = \sum_{n=0}^{\infty} r^n = 1 + r + r^2 + r^3 + \cdots\]
The value of \(r\) determines whether this sum converges (has a finite limit) or diverges.
Understand why this series:
\[\sum_{n=0}^{\infty} r^n\]
converges if and only if:
\[\boxed{|r| < 1}\]
Let’s consider the sum up to the \(N\)-th term:
\[S_N = \sum_{n=0}^{N} r^n = 1 + r + r^2 + \cdots + r^N\]
This has a closed-form expression:
\[S_N = \frac{1 - r^{N+1}}{1 - r} \quad \text{for } r \neq 1\]
To find the sum of the infinite series, take the limit as \(N \to \infty\):
\[S = \lim_{N \to \infty} S_N = \lim_{N \to \infty} \frac{1 - r^{N+1}}{1 - r}\]
If \(|r| < 1\), then:
\[r^{N+1} \to 0 \quad \text{as } N \to \infty\]
So the sum becomes:
\[S = \frac{1}{1 - r}\]
✅ The geometric series converges.
❌ In all cases: the series diverges
\[S = 1 + 0.5 + 0.25 + 0.125 + \cdots = \frac{1}{1 - 0.5} = 2\]
Every term adds less. The sum “flattens out.”
The Z-Transform often gives us geometric series like:
\[\sum_{n=0}^{\infty} (az^{-1})^n\]
This converges only if:
\[|az^{-1}| < 1 \Rightarrow |z| > |a|\]
So, understanding convergence of geometric series = understanding ROC in Z-transforms.
| Condition | Behavior | Result |
|---|---|---|
| \(|r| < 1\) | Terms shrink | Series converges |
| \(|r| \geq 1\) | Terms grow or oscillate | Diverges |
\[\sum_{n=0}^{\infty} r^n = \frac{1}{1 - r} \quad \text{if } |r| < 1\]
For convergence of the geometric series:
\[|az^{-1}| < 1 \Rightarrow |z| > |a|\]
Therefore, the ROC is:
\[\boxed{|z| > |a|}\]
\[x[n] = (0.5)^n u[n]\]
Z-Transform:
\[X(z) = \sum_{n=0}^{\infty} (0.5)^n z^{-n} = \frac{z}{z - 0.5}\]
Region of Convergence:
\[\boxed{|z| > 0.5}\]
| Signal | Z-Transform | ROC |
|---|---|---|
| \(x[n] = a^n u[n]\) | \(\frac{z}{z - a}\) | \(|z| > |a|\) |
| Example: \(a = 0.5\) | \(\frac{z}{z - 0.5}\) | \(|z| > 0.5\) |
Let \(x[n] = a^n u[n]\), where \(|a| < 1\)
\[X(z) = \sum_{n=0}^{\infty} a^n z^{-n} = \frac{1}{1 - az^{-1}}, \quad \text{ROC: } |z| > |a|\]
A difference equation relates input and output values at different time steps.
\[y[n] - a_1 y[n-1] - a_2 y[n-2] = b_0 x[n] + b_1 x[n-1]\]
Common in: - Digital filters (FIR, IIR) - Signal models in ECG, EEG analysis - Implementation in real-time biosignal systems
The Z-Transform turns time shifts into powers of \(z^{-1}\):
| Time Domain | Z-Domain |
|---|---|
| \(x[n]\) | \(X(z)\) |
| \(x[n-k]\) | \(z^{-k} X(z)\) |
| \(y[n-k]\) | \(z^{-k} Y(z)\) |
Given:
\[y[n] - a_1 y[n-1] - a_2 y[n-2] = b_0 x[n] + b_1 x[n-1]\]
Apply \(\mathcal{Z} \{ \cdot \}\):
\[Y(z) - a_1 z^{-1} Y(z) - a_2 z^{-2} Y(z) = b_0 X(z) + b_1 z^{-1} X(z)\]
Group:
\[Y(z)(1 - a_1 z^{-1} - a_2 z^{-2}) = X(z)(b_0 + b_1 z^{-1})\]
Divide both sides:
\[H(z) = \frac{Y(z)}{X(z)} = \frac{b_0 + b_1 z^{-1}}{1 - a_1 z^{-1} - a_2 z^{-2}}\]
Given:
\[y[n] - 0.9 y[n-1] = x[n] - 0.5 x[n-1]\]
Z-Transform:
\[Y(z)(1 - 0.9 z^{-1}) = X(z)(1 - 0.5 z^{-1})\]
Transfer Function:
\[H(z) = \frac{1 - 0.5 z^{-1}}{1 - 0.9 z^{-1}}\]
Let’s analyze \(H(z)\):
Pole-Zero Plot Visualizes system behavior Check for: - Stability (poles inside unit circle) - Frequency shaping
Convert this equation:
\[y[n] = 0.6 y[n-1] + x[n] + x[n-1]\]
Find: - \(H(z)\) - Poles and zeros - Plot them in the Z-plane
Notch Filter
A common issue in biosignal processing is removing power‐line interference (50/60 Hz) from, for example, EEG or ECG signals. A simple digital filter to eliminate 60 Hz interference (assuming a sampling frequency \(f_s = 5000\) Hz) is to place complex‐conjugate zeros at
\[ z = e^{\pm j 2\pi\frac{60}{5000}}. \]
The resulting transfer function can be written as
\[ H(z) = 1 \;-\; 2\cos\!\Bigl(2\pi\frac{60}{5000}\Bigr)\,z^{-1} \;+\; z^{-2}. \]
This \(H(z)\) has zeros at \(e^{\pm j2\pi(60/5000)}\) that precisely cancel the 60 Hz component, thereby implementing a notch filter. Moreover, it is a second‐order FIR filter with symmetric coefficients, which grants it linear phase (important to avoid waveform distortion).
def design_filter(zeros=None, poles=None, gain=1.0):
"""
Diseña un filtro digital a partir de ceros y/o polos y una ganancia.
Parámetros:
- zeros: lista de ceros (raíces del numerador), o None para no incluir
- poles: lista de polos (raíces del denominador), o None para no incluir
- gain: ganancia escalar del filtro
Devuelve:
- b: coeficientes del numerador
- a: coeficientes del denominador
"""
# Si no se pasan ceros, asumimos un FIR trivial (b = [gain])
if zeros:
b = gain * np.poly(zeros)
else:
b = np.array([gain], dtype=float)
# Si no se pasan polos, asumimos sistema FIR (a = [1])
if poles:
a = np.poly(poles)
else:
a = np.array([1.0], dtype=float)
return b, adef simulate_ecg(duration=10.0, fs=500, heart_rate=60):
"""
Simula un ECG sintético basado en la superposición de ondas gaussianas.
Parámetros:
- duration: duración de la señal en segundos
- fs: frecuencia de muestreo en Hz
- heart_rate: frecuencia cardiaca en latidos por minuto
Devuelve:
- t: vector de tiempos
- ecg: señal simulada de ECG en mV
"""
dt = 1 / fs
t = np.arange(0, duration, dt)
rr = 60 / heart_rate # intervalo RR en segundos
# Inicializar señal
ecg = np.zeros_like(t)
# Parámetros de las ondas (posiciones relativas en segundos)
# P wave
p_amp, p_dur, p_delay = 0.25, 0.09, 0.16
# Q wave
q_amp, q_dur, q_delay = -0.05, 0.066, 0.166
# R wave
r_amp, r_dur, r_delay = 1.6, 0.1, 0.166
# S wave
s_amp, s_dur, s_delay = -0.25, 0.066, 0.19
# T wave
t_amp, t_dur, t_delay = 0.35, 0.142, 0.36
# Generar cada latido
for beat_start in np.arange(0, duration, rr):
mask = (t >= beat_start) & (t < beat_start + rr)
tb = t[mask] - beat_start
ecg[mask] += gaussian(tb, p_delay, p_dur / 2, p_amp)
ecg[mask] += gaussian(tb, q_delay, q_dur / 2, q_amp)
ecg[mask] += gaussian(tb, r_delay, r_dur / 2, r_amp)
ecg[mask] += gaussian(tb, s_delay, s_dur / 2, s_amp)
ecg[mask] += gaussian(tb, t_delay, t_dur / 2, t_amp)
return t, ecg+0.1*np.cos(2*np.pi*60*t)# Parámetros de simulación
DURATION = 10 # segundos
FS = 500 # Hz
HR = 70 # latidos por minuto
# Generar señal
t, ecg_signal = simulate_ecg(duration=DURATION, fs=FS, heart_rate=HR)
# Graficar resultado
plt.figure(figsize=(12, 4))
plt.plot(t, ecg_signal, linewidth=1)
plt.title(f'Señal de ECG sintética ({HR} bpm)')
plt.xlabel('Tiempo (s)')
plt.ylabel('Amplitud (mV)')
plt.grid(True)
plt.tight_layout()
plt.show()FIR filters
FIR filters (Finite Impulse Response) are widely used in biomedical processing because they can be designed to have linear phase response, avoiding phase distortion in the filtered signal (which is useful for preserving the morphology of ECG waves, for example).
Filter Design Process
Defining the desired ideal frequency response \(H_d(e^{j\omega})\).
Obtaining the ideal impulse response \(h_d[n]\) as the inverse Fourier transform of \(H_d\).
Truncating \(h_d[n]\) (which is usually infinite or very long) with a window function \(w[n]\) to obtain a realizable FIR filter
\[ h[n] = h_d[n]\,w[n]. \]
You obtain \(h_d[n]\) as the inverse discrete–time Fourier transform of your ideal frequency response \(H_d(e^{j\omega})\). For a low-pass filter with cutoff \(\omega_c\),
\[ H_d(e^{j\omega}) = \begin{cases} 1, & |\omega|\le\omega_c,\\ 0, & \text{otherwise}. \end{cases} \]
By definition of the inverse DTFT,
\[ h_d[n] \;=\; \frac{1}{2\pi}\int_{-\pi}^{\pi} H_d(e^{j\omega})\,e^{j\omega n}\,d\omega \;=\;\frac{1}{2\pi}\int_{-\omega_c}^{\omega_c} e^{j\omega n}\,d\omega. \]
Carry out the integral in two cases:
For \(n\neq 0\):
\[ h_d[n] =\frac{1}{2\pi}\,\frac{e^{j\omega n}}{j\,n}\Biggr|_{-\omega_c}^{\omega_c} =\frac{1}{2\pi}\,\frac{e^{j\omega_c n}-e^{-j\omega_c n}}{j\,n} =\frac{\sin(\omega_c\,n)}{\pi\,n}. \]
For \(n=0\):
\[ h_d[0] =\frac{1}{2\pi}\int_{-\omega_c}^{\omega_c} 1\,d\omega =\frac{2\,\omega_c}{2\pi} =\frac{\omega_c}{\pi}. \]
Putting both together,
\[ h_d[n] =\begin{cases} \dfrac{\sin(\omega_c\,n)}{\pi\,n}, & n\neq 0,\\[1em] \dfrac{\omega_c}{\pi}, & n=0. \end{cases} \]
If your cutoff is specified in Hz, \(f_c\), and sampling rate is \(f_s\), then
\[ \omega_c = 2\pi\,\frac{f_c}{f_s}, \]
so you can write
\[ h_d[n] =\frac{\sin\!\bigl(2\pi \frac{f_c}{f_s}\,n\bigr)}{\pi\,n} =\;2\;\frac{f_c}{f_s}\;\frac{\sin\!\bigl(2\pi \frac{f_c}{f_s}\,n\bigr)}{2\pi \frac{f_c}{f_s}\,n} =2\frac{f_c}{f_s}\,\mathrm{sinc}\!\Bigl(2\frac{f_c}{f_s}\,n\Bigr), \]
where we define the normalized sinc as \(\mathrm{sinc}(x)=\frac{\sin(\pi x)}{\pi x}\).
The typical characteristics of classic windows are:
Rectangular: narrowest main lobe (width ≈ 4π/M rad) but highest sidelobes (first sidelobe ≈ –13 dB, stop-band attenuation ≈ 21 dB). It gives the steepest transition for a given filter order, at the expense of poorer stop-band rejection.
Bartlett (triangular): somewhat wider main lobe (≈ 8π/M), sidelobes ≈ –25 dB.
Hann: main lobe ≈ 8π/M, sidelobes ≈ –31 dB (better rejection than rectangular but smoother transitions).
Hamming: main lobe ≈ 8π/M, sidelobes ≈ –41 dB (minimum stop-band attenuation ≈ 53 dB). Very popular for its good compromise between transition width and stop-band attenuation.
Blackman: wider main lobe (≈ 12π/M) but very low sidelobes (≈ –57 dB, attenuation ≈ 74 dB).
Kaiser: allows selection of a parameter β to control sidelobe attenuation, offering flexibility. Approximately, to achieve A dB of attenuation,
\[ \beta \approx 0.1102\,(A - 8.7)\quad(\text{for }A>50), \]
and the normalized transition width Δω relates to the order M and β by
\[ M \approx \frac{A - 8}{2.285\,\Delta\omega} \]
These formulas stem from Kaiser’s approximations and help in sizing the filter.
(-100.0, 5.0)
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal.windows import (
boxcar, # Rectangular
bartlett, # Triangular
hann, # Hann
hamming, # Hamming
blackman, # Blackman
kaiser # Kaiser
)
# Parameters
M = 64 # window length
beta = 8.6 # Kaiser parameter for moderate sidelobe attenuation
nfft = 512 # FFT length for frequency response
fs = 1.0 # normalized sampling rate
# Generate windows
windows = {
'Rectangular': boxcar(M),
'Bartlett': bartlett(M),
'Hann': hann(M),
'Hamming': hamming(M),
'Blackman': blackman(M),
f'Kaiser (β={beta})': kaiser(M, beta)
}
# Time-domain plot
plt.figure(figsize=(8, 4))
for name, w in windows.items():
plt.plot(np.arange(M), w, label=name)
plt.title('Window Functions — Time Domain')
plt.xlabel('Sample index n')
plt.ylabel('Amplitude')
plt.legend(loc='upper right')
plt.grid(True)
plt.tight_layout()
# Frequency-domain plot
plt.figure(figsize=(8, 4))
for name, w in windows.items():
# Compute normalized frequency response
W = np.fft.fft(w, nfft)
W_mag = 20 * np.log10(np.abs(W) / np.max(np.abs(W)))
freqs = np.fft.fftfreq(nfft, d=1/fs)
# Only plot 0 ≤ f ≤ 0.5 (normalized Nyquist)
idx = np.logical_and(freqs >= 0, freqs <= 0.5)
plt.plot(freqs[idx], W_mag[idx], label=name)
plt.title('Window Functions — Frequency Response')
plt.xlabel('Normalized Frequency (cycles/sample)')
plt.ylabel('Magnitude (dB)')
plt.ylim(-100, 5)
plt.legend(loc='upper right')
plt.grid(True)
plt.tight_layout()
plt.show()N = 101
fc = 30
# 2. Compute the ideal impulse response h_d[n] of a low-pass filter
n = np.arange(N)
M = (N - 1) / 2
# Using the normalized sinc: sinc(x) = sin(pi*x)/(pi*x)
h_d = (2 * fc / FS) * np.sinc(2 * fc * (n - M) / FS)
# 3. Choose a window w[n] (here: Hamming) and truncate
w = np.hamming(N)
h = h_d * w
# 4. Normalize to ensure unity gain at DC
h /= np.sum(h)
# 6. Filter it (only allowed library call)
y = sig.lfilter(h, 1.0, ecg_signal)
# 7. Plot input vs. output
plt.figure(figsize=(8, 4))
plt.plot(t, ecg_signal, label='Original signal')
plt.plot(t, y, label='Filtered signal', linewidth=2)
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Low-pass FIR (window method) – Cutoff 100 Hz')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()Designing an analog filter translates desired signal filtering characteristics into a physical electronic circuit.
Key Stages:
We’ll illustrate steps 1, 3, and aspects of 4/5 using Python’s scipy.signal.
This is the most crucial first step.
gpass)gstop)Let’s define specs for an analog low-pass filter. Frequencies for analog filters in scipy are typically in rad/s.
Passband edge: 6283.19 rad/s (1000 Hz)
Stopband edge: 9424.78 rad/s (1500 Hz)
Max Passband Loss: 1.0 dB
Min Stopband Attenuation: 40.0 dB
Select based on trade-offs (roll-off, ripple, phase).
For our Python examples, we’ll use Butterworth.
The Butterworth filter is maximally flat in the passband.
The normalized prototype has a cutoff frequency \(\omega_c = 1\) rad/s.
Poles are placed on the unit circle in the \(s\)-plane at angles: \[s_k = e^{j\pi\frac{2k + n - 1}{2n}}, \quad k = 1,2,\ldots,n.\]
The Butterworth polynomial of order \(n\) is: \[B_n(s) = \prod_{k=1}^{n} (s - s_k).\]
The normalized \(n\)th-order lowpass prototype is: \[H_{LP}(s) = \frac{1}{B_n(s)} = \frac{1}{\prod_{k=1}^n (s - s_k)}.\] This yields a cutoff at \(|H_{LP}(j1)| = 1/\sqrt{2}\) (the –3 dB point).
Apply the frequency transformation \(s \to 1/s\) and multiply by \(s^n\): \[H_{HP}(s) = H_{LP}(1/s)\,s^n = \frac{s^n}{B_n(1/s)} = \frac{s^n}{\prod_{k=1}^n (1/s - s_k)} = \frac{\prod_{k=1}^n(-s_k)\;s^n}{\prod_{k=1}^n(s - s_k^{-1})}.\] This prototype has a normalized cutoff at \(\omega=1\) rad/s in the highpass sense.
Use the bandpass transformation: \[s \to \frac{s^2 + \Omega_0^2}{B\,s},\] where \(\Omega_0 = \sqrt{\omega_1\omega_2}\) and \(B = \omega_2 - \omega_1\).
For a normalized prototype, we set \(\Omega_0 = 1\), \(B = 1\), giving: \[H_{BP}(s) = H_{LP}\left(\frac{s^2 + 1}{s}\right) = \frac{(s)^n}{\prod_{k=1}^n\left(\frac{s^2 + 1}{s} - s_k\right)} = \frac{s^n}{\prod_{k=1}^n\left(s^2 + 1 - s\,s_k\right)}.\]
Use the bandstop transformation: \[s \to \frac{s}{\frac{s^2 + 1}{s}} = \frac{s^2}{s^2 + 1},\] or equivalently \(s \to (B\,s)/(s^2 + \Omega_0^2)\) with \(B=1\), \(\Omega_0=1\). Then: \[H_{BS}(s) = H_{LP}\left(\frac{s}{s^2 + 1}\right) = \frac{1}{B_n\left(\frac{s}{s^2+1}\right)} = \frac{1}{\prod_{k=1}^n\left(\frac{s}{s^2+1} - s_k\right)} = \frac{(s^2+1)^n}{\prod_{k=1}^n\left(s - s_k(s^2+1)\right)} = \frac{(s^2+1)^n}{\prod_{k=1}^n\left(s - s_ks^2 - s_k\right)}.\]
Estimating the filter order (\(N\)) is a crucial first step in analog filter design after defining specifications. The order determines the filter’s complexity and steepness of its frequency response roll-off.
These equations are for Butterworth filters. The calculated \(N\) is a real number and must be rounded up to the next integer.
Common Variables: * \(N\): Filter order (integer, result of formula is rounded up). * \(A_{min}\): Minimum stopband attenuation (dB, positive value, e.g., 40 dB). * \(A_{max}\): Maximum passband attenuation/ripple (dB, positive value, e.g., 1 dB). * \(\log_{10}\): Base-10 logarithm. * All \(\Omega\) values are angular frequencies in rad/s.
\[N \ge \frac{\log_{10}\left(\frac{10^{A_{min}/10}-1}{10^{A_{max}/10}-1}\right)}{2 \log_{10}\left(\frac{\Omega_s}{\Omega_p}\right)}\]
Explanation of LP-Specific Variables: * \(\Omega_p\): Passband edge angular frequency. The frequency up to which signals pass with at most \(A_{max}\) attenuation. * \(\Omega_s\): Stopband edge angular frequency. The frequency from which signals are attenuated by at least \(A_{min}\). * For a low-pass filter, \(\Omega_s > \Omega_p\). The term \(\frac{\Omega_s}{\Omega_p}\) is the transition ratio.
\[N \ge \frac{\log_{10}\left(\frac{10^{A_{min}/10}-1}{10^{A_{max}/10}-1}\right)}{2 \log_{10}\left(\frac{\Omega_p}{\Omega_s}\right)}\]
Explanation of HP-Specific Variables: * \(\Omega_p\): Passband edge angular frequency. The frequency from which signals pass with at most \(A_{max}\) attenuation. * \(\Omega_s\): Stopband edge angular frequency. The frequency down to which signals are attenuated by at least \(A_{min}\). * For a high-pass filter, \(\Omega_p > \Omega_s\). The term \(\frac{\Omega_p}{\Omega_s}\) is the transition ratio (inverted compared to LP to keep it >1 for the logarithm).
The order \(N_{BP}\) is the same as an equivalent Low-Pass Prototype. \[N_{BP} \ge \frac{\log_{10}\left(\frac{10^{A_{min}/10}-1}{10^{A_{max}/10}-1}\right)}{2 \log_{10}\left(\Omega_{s,LP_equiv}\right)}\]
Explanation of BP-Specific Variables: * \(\Omega_{p1}, \Omega_{p2}\): Passband edge angular frequencies (\(\Omega_{p1} < \Omega_{p2}\)). * \(\Omega_{s1}, \Omega_{s2}\): Stopband edge angular frequencies (\(\Omega_{s1} < \Omega_{p1}\) and \(\Omega_{s2} > \Omega_{p2}\)). * \(\Omega_0 = \sqrt{\Omega_{p1}\Omega_{p2}}\): Geometric center frequency of the passband. * \(B = \Omega_{p2} - \Omega_{p1}\): Bandwidth of the passband. * \(\Omega_{s,LP_equiv} = \min\left( \left| \frac{\Omega_{s1}^2 - \Omega_0^2}{\Omega_{s1} B} \right|, \left| \frac{\Omega_{s2}^2 - \Omega_0^2}{\Omega_{s2} B} \right| \right)\) * This \(\Omega_{s,LP_equiv}\) is the stopband edge of an equivalent low-pass prototype whose passband edge is normalized to 1 rad/s. It represents the more stringent of the two possible transitions from the passband edges to the stopband edges.
The order \(N_{BS}\) is the same as an equivalent Low-Pass Prototype. \[N_{BS} \ge \frac{\log_{10}\left(\frac{10^{A_{min}/10}-1}{10^{A_{max}/10}-1}\right)}{2 \log_{10}\left(\Omega_{trans,LP_equiv}\right)}\]
Explanation of BS-Specific Variables: * \(\Omega_{s1}, \Omega_{s2}\): Stopband edge angular frequencies defining the rejected band (\(\Omega_{s1} < \Omega_{s2}\)). \(A_{min}\) is the attenuation in this band. * \(\Omega_{p1}, \Omega_{p2}\): Passband edge angular frequencies outside the stopband (\(\Omega_{p1} < \Omega_{s1}\) and \(\Omega_{p2} > \Omega_{s2}\)). \(A_{max}\) is the ripple in these passbands. * \(\Omega_0 = \sqrt{\Omega_{s1}\Omega_{s2}}\): Geometric center frequency of the stopband. * \(B = \Omega_{s2} - \Omega_{s1}\): Bandwidth of the stopband. * \(\Omega_{p,LP_i} = \left| \frac{\Omega_{pi} B}{\Omega_{pi}^2 - \Omega_0^2} \right|\) for \(i=1\) (using \(\Omega_{p1}\)) and \(i=2\) (using \(\Omega_{p2}\)). * These are the passband edges of an equivalent low-pass prototype whose stopband edge is normalized to 1 rad/s. * \(\Omega_{trans,LP_equiv} = \frac{1}{\min(\Omega_{p,LP1}, \Omega_{p,LP2})}\) * This is the effective transition ratio (stopband edge / passband edge) for the equivalent low-pass prototype.
Calculated based on specs and chosen filter type. scipy.signal.buttord (and similar for other types) can find this for analog filters.
Required Butterworth Filter Order (N): 14
Butterworth Natural Frequency (Wn_buttord): 6593.83 rad/s
Wn_buttord is the natural angular frequency (\(\Omega_n\)) for the filter. For Butterworth, this is the -3dB cutoff frequency \(\Omega_c\).Conceptually: Find normalized LP prototype \(\rightarrow\) Denormalize \(\rightarrow\) Transform.
scipy.signal functions like butter, cheby1, etc., (with analog=True) perform these steps internally to give the final analog filter transfer function \(H_a(s)\).
\(H_a(s)\) can be represented as: 1. Numerator (\(b\)) and Denominator (\(a\)) coefficients: \(H_a(s) = \frac{b_0 s^M + b_1 s^{M-1} + \dots + b_M}{a_0 s^N + a_1 s^{N-1} + \dots + a_N}\) 2. Zeros (\(z\)), Poles (\(p\)), and Gain (\(k\)): \(H_a(s) = k \frac{\prod (s-z_i)}{\prod (s-p_j)}\)
Using signal.butter with analog=True. The Wn parameter here should be the cutoff frequency \(\Omega_c\). For buttord, the returned Wn_buttord is exactly this \(\Omega_c\).
Analog Filter Coefficients (Numerator b, Denominator a):
b_analog = [2.93718174e+53]
a_analog = [1.00000000e+00 5.88921844e+04 1.73414469e+09 3.37531810e+13
4.84169654e+17 5.40639098e+21 4.84125527e+25 3.52958830e+29
2.10491137e+33 1.02201934e+37 3.97946840e+40 1.20619642e+44
2.69441500e+47 3.97843852e+50 2.93718174e+53]
Analog Filter Zeros, Poles, Gain (ZPK):
Zeros (z_analog) = []
Poles (p_analog) = [ -738.27500968+6552.37193897j -2177.8048373 +6223.80864963j
-3508.13043664+5583.15760623j -4662.54372722+4662.54372722j
-5583.15760623+3508.13043664j -6223.80864963+2177.8048373j
-6552.37193897 +738.27500968j -6552.37193897 -738.27500968j
-6223.80864963-2177.8048373j -5583.15760623-3508.13043664j
-4662.54372722-4662.54372722j -3508.13043664-5583.15760623j
-2177.8048373 -6223.80864963j -738.27500968-6552.37193897j]
Gain (k_analog) = 293718173729774468825490597908032846661500110135885824.00
Use scipy.signal.freqs to compute the frequency response of an analog filter given its \(b,a\) coefficients.
(100.0, 7500)
(-60.0, 5.0)
(100.0, 7500)
Convert the mathematical \(H_a(s)\) into a physical electronic circuit. This step is not done by scipy.signal. It requires circuit design knowledge.
This is also outside the scope of scipy.signal but critical for hardware implementation.
scipy.signal: Excellent for determining filter parameters (\(N, \Omega_c\)) and transfer functions (\(b,a\) or \(z,p,k\)), and plotting responses.Analog filter design is a multi-step process:
scipy.signal in Python is a powerful tool for the mathematical parts: determining order, calculating transfer function coefficients/poles/zeros, and analyzing frequency response for various analog filter prototypes.gpass (dB)gstop (dB)Sampling Frequency: 10000 Hz
Normalized Passband Edge: 0.200 (pi rad/sample)
Normalized Stopband Edge: 0.300 (pi rad/sample)
Max Passband Loss: 1 dB
Min Stopband Attenuation: 40 dB
scipy.signal functions for IIR design.scipy.signal.buttord to find the minimum order and cutoff frequency for a digital Butterworth filter meeting the specs.
Minimum Filter Order (N): 12
Butterworth Natural Frequency (Wn): 0.211 (pi rad/sample)
analog=False specifies we want parameters for a digital filter.Wn is the digital cutoff frequency (scalar for LP/HP, tuple for BP/BS) required by the butter function. It’s often close to wp.scipy.signal.butter (or cheby1, cheby2, ellip) to get the filter coefficients \(b\) (numerator) and \(a\) (denominator) of \(H(z)\).
Filter Coefficients:
Numerator (b): [0.0e+00 0.0e+00 1.0e-05 4.0e-05 1.0e-04 1.6e-04 1.9e-04 1.6e-04 1.0e-04
4.0e-05 1.0e-05 0.0e+00 0.0e+00]
Denominator (a): [ 1.000000e+00 -6.930840e+00 2.271843e+01 -4.632962e+01 6.523097e+01
-6.662310e+01 5.050575e+01 -2.858485e+01 1.197050e+01 -3.612930e+00
7.452400e-01 -9.424000e-02 5.520000e-03]
output='ba' gives numerator/denominator coefficients. Other options: 'sos' (second-order sections - numerically better), 'zpk' (zeros, poles, gain).scipy.signal.freqz.(0.0, 5000.0)
(-60.0, 5.0)
(0.0, 5000.0)
b, a define the difference equation.scipy.signal.lfilter to apply the filter to a signal.(-1.5, 1.5)
buttord (digital).butter (digital).freqz.lfilter (uses the difference equation implicitly).